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An atomistic method of calculating the spin-lattice relaxation times (Ti) is presented for donors 
in silicon nanostructures comprising of millions of atoms. The method takes into account the full 
band structure of silicon including the spin-orbit interaction. The electron-phonon Hamiltonian, 
and hence the deformation potential, is directly evaluated from the strain-dependent tight-binding 
Hamiltonian. The technique is applied to single donors and donor clusters in silicon, and explains 
the variation of Ti with the number of donors and electrons, as well as donor locations. Without any 
adjustable parameters, the relaxation rates in a magnetic field for both systems are found to vary 
as B 5 in excellent quantitative agreement with experimental measurements. The results also show 
that by engineering electronic wavefunctions in nanostructures, T\ times can be varied by orders of 
magnitude. 

PACS numbers: 71.55.Cn, 03.67.Lx, 85.35.Gv, 71.70.Ej 


Due to the extremely long spin coherence times, in 
some cases exceeding seconds mm, and the existing in¬ 
dustrial fabrication infrastructure, silicon is well-suited to 
be an outstanding platform for semiconductor quantum 
computer technology nu. Qubits hosted by donors in 
silicon [3] have some added advantages as they are read¬ 
ily available few-electron systems with a rich electronic 
structure and can form identical qubits [9]. In the last 
few years, several key experimental milestones have been 
achieved in dopant based quantum computing, including 
the demonstration of electron m and nuclear m spin 
qubits, single spin read-out and initialization m he 
and the observation of spin blockade and exchange to¬ 
wards two qubit coupling Hang. Recent advances in 
Scanning Tunneling Microscope (STM) lithography has 
enabled placement of single donors with atomic scale pre¬ 
cision m, with the result that various functional units 
such as quantum wires m, single electron transistors 
(SET) TS! (18|. and quantum dots |19j can all be realized 
in-plane with densely packed donor islands. The STM 
approach provides the fabrication precision needed to de¬ 
velop test-bed quantum chips for the demonstration of 
quantum algorithms in a solid-state quantum computer. 

One of the two most important timescales for a spin 
qubit is the spin-lattice relaxation time (Ti). Recent ex¬ 
periments have measured T\ times in a single donor and 
in a few-donor cluster indicating shorter T\ times in the 
latter nans. Previous theoretical works exist in the lit¬ 
erature qualitatively describing two different spin relax¬ 
ation mechanisms in a bulk donor system IIBEI]. How¬ 
ever, a comprehensive quantitative theory which com¬ 
bines all the different mechanisms under a unified frame¬ 
work and accounts for the local inhomogeneous environ¬ 
ment of the donors in a realistic nanostructure is still 


lacking. Moreover, there is no theoretical work yet to ex¬ 
plain the measured T\ times in densely packed donor clus¬ 
ters. In this letter, we present a comprehensive approach 
to compute the T\ times in single donors and donor clus¬ 
ters in silicon nanostructures based on the self-consistent 
atomistic tight-binding (TB) method. The computed T\ 
times can explain the experimental results of Refs Hana 
without any adjustable parameters. The T\ times are 
found to depend strongly on the size and shape of the 
electronic wavefunctions, which suggests that quantum 
confinement plays an important role in the relaxation 
process. The calculations also provide an insight into 
how the T\ times can be engineered by several orders of 
magnitude in these quantum devices. 



FIG. 1: a) Schematic plot of a 4P donor cluster in silicon. 
The cluster hosts an up spin electron with energy E n which 
relaxes to the down spin state E n > by emitting a phonon with 
energy huj q . b) An STM template of a 4P cluster showing 
example dimer locations 13]. c) The computed self-consistent 
probability density of the outermost electron in a 4P donor 
cluster with 5 electrons (4P5e). 

Single shot measurements of the donor-bound elec- 
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tron spin performed in ion-implanted metal-oxide- 
semiconductor (MOS) devices yielded a T\ time of 2.3 
s at B = 2 T and a B 5 dependency of the relaxation 
rate (1/Ti) [12 . Similar T\ measurements have since 
been performed in a STM patterned donor cluster in an 
all epitaxially engineered device with atomic scale preci¬ 
sion, giving T\ time of 0.4 s at B = 2 T and a similar 
B 5 dependency of the relaxation rate [13]. The donor 
clusters are composed of several donors with an average 
separation less than or equal to the single donor Bohr 
radii (Fig. lb) which can host several electrons. Such 
a cluster provides for additional addressability within an 
array of qubits m- In general, the T\ times of donors 
in realistic devices and nanostructures are likely to be 
influenced by the inhomogeneous environment, and need 
to be understood from an atomic scale theory. 

The full TB Hamiltonian of the silicon and the P 
atoms was represented with a 20-orbital sp 3 d 5 «s* basis 
per atom including nearest-neighbor and spin-orbit inter¬ 
actions |23ll25l . A donor was represented by a Coulomb 
potential of a positive charge screened by the dielectric 
constant of Si and subjected to an onsite cut off potential 
[26] [27]. The model reproduces the full energy spectrum 
of a single donor including valley-orbit splitting [22] , and 
also captures the single donor hyperfine [26] and spin- 
orbit Stark effects [28] in close agreement with experi¬ 
ments. Recent STM imaging of the donor wavefunction 
also shows excellent agreement with the tight-binding 
wavefunction at the atomic scale m- The magnetic field 
is represented by a vector potential in a symmetric gauge 
and entered through a Peierls substitution. To capture 
multi-electron occupation of the donor clusters, a self- 
consistent Hartree method was employed, in which the 
electronic charge density was computed from the N low¬ 
est energy occupied wavefunctions, n(r) = |4q(r)| 2 . 

Solving the potential due to n(r) self-consistently with 
the tight-binding Hamiltonian until convergence enables 
us to obtain the binding energies and the wavefunctions 
of the few electrons bound to the donor cluster. This 
method has also successfully reproduced the experimen¬ 
tal D~ (i.e., the 2nd electron bound to a single P) binding 
energy USES]. Exchange and correlation terms based 
on local density approximation are added to the Hartree 
potential m- The same methodology has been used to 
reproduce experimentally measured addition energies of 
multi-donor clusters m The TB Hamiltonian of about 
1.4 million atoms including the Hartree potential is then 
solved by a parallel Block Lanczos algorithm to obtain 
the relevant lowest energy wavefunctions. 

For the relaxation times computed in this work, we 
assume that an electron has been loaded to the ground 
state of a single donor or a donor cluster. In case of a bulk 
P donor, this represents the A\ state at -45.6 meV be¬ 
low the conduction band [32] . This binding energy can 
vary with donor depths and fields in a realistic device 
[22]. In addition, the binding energy in a donor cluster 


can be sensitive to the donor numbers, electron numbers, 
and donor locations. Experimentally, the electron can be 
loaded in the ground state of the donor/cluster by bring¬ 
ing the Fermi level of an electron reservoir in resonance 
with the Zeeman split ground state, as demonstrated in 
Refs nans]. 

The relaxation rate 1/Ti, for an electron-phonon in¬ 
teraction Hamiltonian H ep can be obtained by Fermi’s 
Golden Rule, 

1 27r 

= -r\(n',n q + l\H ep \n,n q )\ 2 5(E n - E n > - hu q ) ( 1 ) 
1 1 ri 

where n and n' are the up and down spin electronic states 
with energy E n and E n > respectively, uo q is the angular 
frequency of the emitted phonon, n q and n q + 1 are the 
initial and final phonon states with wavevector q. In 
the B-field range of experimental interest [lJJ [13] . the 
Zeeman splitting is less than 1 meV. As a result, only 
acoustic phonons contribute to the spin relaxation pro¬ 
cess. H ep then depends on the deformation potential of 
the crystal (i, j representing each of the three Carte¬ 
sian directions) and the strain tensor components [33] 
(both of which are position dependent in the atomistic 
TB method), and is given by 

Hep = Ztjj Uij (2) 

hj 

To evaluate S^*, we use the relation and 

compute the total change in the electron-phonon Hamil¬ 
tonian A H ep due to an infinitesimal uniform lattice strain 
represented by A Uij = (a small arbitrary constant). 
Since A H ep can be expressed as a change in the electronic 
TB Hamiltonian H e under a crystal deformation caused 
by [33} [34] , is given as 

1 - J y = {He{uij) — Hem/Uij (3) 

The strain dependent TB Hamiltonian H e ( Uij ) expresses 
the TB matrix elements as functions of inter-atomic bond 
lengths and distortion angles depending on the relative 
positions of pairs of atoms in the lattice [35] . This 
method of incorporating local strain in the TB Hamil¬ 
tonian is well-established and has been shown to repro¬ 
duce various experimental results [23]. Although we have 
considered all 6 components of , we have found the off- 
diagonal terms to be small for single donor states located 
near the conduction band. 

Furthermore, Uij can also be expressed in terms of 
phonon creation and annhilation operators, and a q 
respectively, as 

Uij = \ F)( 2 A Mw + eqj<7i)(a+ exp[i(g • r)] + 

q ^ q 

a q exp[-i(q ■ r)]) 

( 4 ) 
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where V is the volume of the crystal, p the mass density, 
and e q the phonon polarization unit vector. Using eq. 2 
and 4, the matrix element of H ep can be expressed as, 

(n',n q + l\H ep \n,n q ) = 4 vS + 1 

_ 9 9 . (5) 

2^*( e qi<?j + e qjQi)( n '\ exp[i(q • r)]Sij|n) 

i,j 

In eq. 5, we have used p = 2330 kg/m 3 , while uj q is 
obtained from the electron Zeeman energy E z = fuj q . 
The phonon number n q satisfies the Boltzmann distri¬ 
bution. We choose a temperature range below lOOmK, 
which guarantees us to be in the low-temperature regime 
where (n q + 1) ^ 1 [12]. Above IK, (n q + 1) ~ n g , and 
1/Ti varies as B A [20] ‘21 . The polarization vector e q 
takes into account the three phonon modes: one longi¬ 
tudinal and two transverse. Since the Zeeman splitting 
energy is very small (< 1 meV), linear and isotropic bulk 
Si phonon dispersion is assumed. The phonon wavevector 
q is then evaluated as uj q = v s q , where v s is the speed of 
sound in Si, and taken to be 8480 m/s for the longitudi¬ 
nal mode and 5860 m/s for transverse modes. The same 
constants were also used to interpret the experimental 
data dang. While local vibrational modes have been 
observed for P atoms in Si [36] . the energy corresponds 
to the mode frequency is at least an order of magnitude 
larger than the energy range in interest. The measured 
T\ values do not deviate from the B 5 behavior indicates 
that the local vibrational modes do not contribute signif¬ 
icantly to the spin-relaxation process. 

Fig. 2 shows the spin-lattice relaxation rates, 1/Ti, of 
a single P donor and a 4P donor cluster as a function 
of B-field. The red squares are the measured rates for a 
single donor from Ref. m, whereas the blue triangles 
are the measured rates for a STM patterned few donor 
cluster from Ref. m- Experimentally, the exact number 
of donors in the cluster was unknown, but estimated to be 
between 2 and 5 based on STM images. From transport 
measurements, it was also expected that during the spin 
readout step at least three electrons occupied the donor 
cluster, while in total seven charge transitions on the 
donor cluster were observed [13 . 

The red solid line in Fig. 2 represents the calculations 
performed in this work for a bulk P donor. The calcu¬ 
lated rates show a B 5 dependence of 1/Ti, and also yield 
similar magnitudes of T\ (~2.5 s at B = 2 T) as the ex¬ 
periment. The B field in this calculation is applied along 
the [110] direction, consistent with the experiment |12(. 
To understand the effect of donor number and electron 
number on T\ , we have simulated donor clusters compris¬ 
ing of 2 to 4 donors with various electron occupation. In 
Fig. 2, we show the results of the 4P cluster with 1, 3 
and 5 bound electrons (the green, cyan, and the blue solid 
lines, respectively). While the B 5 dependency holds in 
all cases, the rates vary considerably in magnitude, and 


increase with the number of electrons. Our calculations 
show that higher measured relaxation rates of the cluster 
come from a 5e occupation in a 4P cluster, which is also 
consistent with the experimental finding of the electron 
number being > 3 [15] . 

We have intentionally chosen an odd number of elec¬ 
trons because the relaxation between a net 1/2 and -1/2 
spin is assumed, which requires an unpaired number of 
electrons. This is also consistent with the experimen¬ 
tal measurements, where no spin read-out signal was ob¬ 
served for alternate electron occupation in the cluster 
m ■ The calculations also reveal the startling fact that 
if we have only one electron in a 4P cluster, the relax¬ 
ation rate is actually smaller than the bulk P, and the T\ 
times can be increased from few seconds to hundreds of 
seconds. The physical reason that determines the mag¬ 
nitude of the T\ times is also investigated in this work. 



FIG. 2: Spin-lattice relaxation rates of a single P donor and 
a 4P donor cluster as a function of B-field. The red squares 
and blue triangles show the measured data for a single donor 
[12] and a few-donor cluster m, respectively. The solid lines 
show the TB calculation results for lPle (red), 4Ple (green), 
4P3e (cyan) and 4P5e (blue). 1/Ti varies as B 5 for all cases. 
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FIG. 3: The computed probability density (in log scale) of 
the outermost electron in various donor clusters. The plots 
show a 2D cut through the center of the clusters. 

To understand the impact of electron number on the 
T\ times observed we determined the electron probabil- 
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ity densities of the outermost electrons for varying clus¬ 
ter sizes and electron number in Fig. 3. The size and 
shape of the wavefunctions depends on the number of 
donors and electrons and their locations within the clus¬ 
ter. For the same number of electrons, more donors re¬ 
sult in a more tightly bound wavefunction because of 
the stronger potential of the larger number of positively 
charged donor cores. For the same number of donors, as 
more electrons are added, the wavefunction spreads out 
more as the donor core becomes more strongly screened. 
Electron-electron repulsion also causes the wavefunctions 
of the outermost electron to spread out more. We have 
extracted the Bohr radii of these wavefunctions by fitting 
an exponential decay function to the tail of probability 
densities (|T(r)| 2 ) along the x-axis through the center of 
the clusters. 

Fig. 4a shows the relaxation rates (at B = 2 T) as 
a function of the Bohr radii for the same clusters as in 
Fig. 3. It is observed that donor clusters with larger 
Bohr radii (i.e., those clusters with more electrons and 
fewer donors) result in higher relaxation rates. Since the 
acoustic phonon wavelength corresponding to a Zeeman 
energy of 0.2 meV (at B = 2 T) is about 100 nm, the 
phonon wavelength is much larger than the electronic 
wavelengths in this system. A larger electronic wave- 
function therefore interacts with the phonons more |35]. 
Perhaps more importantly, a larger wavefunction also im¬ 
plies less quantum confinement in the system, and hence 
a smaller energy gap between excited and ground state 
(i.e. the valley-orbit gap) |39|, which relates well to the 
Hasegawa theory for a bulk donor [20] . in which T\ in¬ 
creases with this energy gap. Since our calculations show 
a valley-orbit gap ranging from 30 meV to 5 meV as clus¬ 
ter wavefunction increases in radii, we expect a strong 
dependence of T\ on radii as well. All the calculations 
presented here also include exchange and correlation ef¬ 
fects. We observed that the inclusion of the exchange 
and correlation effects results in slightly larger wavefunc¬ 
tions, as the electrons experience greater net repulsion. 
This causes the relaxation rates to increase slightly and 
move closer to the experimental data in Fig. 2. 

Within a lithographic template for a donor cluster, 
there can be some uncertainties in the exact locations 
of the donor [16.. However, if the cluster only comprises 
of a few donors, all the positional configurations can be 
enumerated, and the most compact and the most dis¬ 
persed donor clusters can be identified. Since it is com¬ 
putationally time intensive to simulate all possible donor 
configurations within a cluster, we have simulated the 
two extreme cases for donor clusters of 1 to 4 donors (Fig. 
4c). In Fig. 4b, we have plotted the relaxation rates as 
a function of electron number for different donor clus¬ 
ters with these two positional configurations (marked by 
crosses within an ellipse). As expected, the T\ times have 
some associated variations with donor locations within a 
cluster. However, the dependency on the total donor and 



FIG. 4: (a) The relaxation rates of various donor clusters 
at B = 2 T as a function of the computed Bohr radii of 
the wavefunction. (b) The relaxation rates of various donor 
clusters at B = 2 T with different number of electrons. The 
different points within an ellipse represent the variations in 
Ti obtained with different donor locations within the cluster, 
(c) Exact donor locations for the donor clusters in (b). 


total electron numbers is stronger. This suggests that Ti 
measurements can be used to infer information about 
donor and electron numbers in STM patterned donor de¬ 
vices, as a non-invasive metrology technique. Such in¬ 
formation can be useful for engineering pulses to control 
single- and two-qubit operations in experiments. 

Previous effective mass calculations of Hasegawa [20] 
and Roth [21 predicted two different spin relaxation 
mechanisms due to an effective g-factor shift. Hasegawa’s 
mechanism predicts this effective g-factor shift due to the 
strain induced redistribution of the donor wavefunction 
among the six conduction band valleys [20] . while Roth’s 
mechanism predicts a single valley g-factor shift due to 
a strain induced mixing of higher conduction bands [151: . 
Both mechanisms inherently depend on the spin-orbit in¬ 
teraction in silicon, which reduces the bulk donor g-factor 
slightly below 2. Both Hasegawa and Roth’s theory show 
a B 5 dependency of 1/Ti for a bulk donor, but are rather 
qualitative in nature. The above approaches are also 
limited to a single bulk donor for which Kohn-Luttinger 
wavefunctions can be used |32|. Our TB method captures 
both mechanisms under the same framework, as a full 
bandstructure description is used from the atomic orbital 
basis including spin-orbit interaction. The method is also 
general and applies to any nanostructures and semicon¬ 
ductors for which accurate TB models can be developed. 

In conclusion, we have presented an atomistic approach 
to calculate the phonon induced spin lattice relaxation 
rates in donors in silicon. The T\ times agree very well 
with recently measured values on single donors and donor 
clusters, and help to explain the variation of T\ with the 
numbers of donors and electrons, and the donor locations. 
The values of T\ were found to have a strong dependency 
on the size of the electronic wavefunctions. This also pro¬ 
vides a way to engineer larger T\ times by using donor 
clusters with large number of P donors and single elec- 
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tron. An atomistic description of the T\ times and their 
variations in inhomogeneous environment provides cru¬ 
cial information in the design of silicon qubits. 
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